#include <iostream>
#include <ale.hpp>


//template<int Dim,int Order> const uint16_type AdvReacDiff<Dim>::Order;

template<int Dim,int Order>
typename AdvReacDiff<Dim,Order>::mesh_ptrtype
AdvReacDiff<Dim,Order>::createMesh( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    GmshTensorizedDomain<convex_type::nDim,
                         convex_type::nOrder,
                         convex_type::nDim,
                         Simplex> td;
    td.setCharacteristicLength( meshSize );
    td.setX( std::make_pair( 0, 5 ) );
    td.setY( std::make_pair( -1, 1 ) );
    std::string fname = td.generate( convex_type::name().c_str() );

    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );


    return mesh;
} // AdvReacDiff::createMesh

#include <life/lifediscr/ale.hpp>
template<int Dim,int Order>
void
AdvReacDiff<Dim,Order>::run()
{
    if ( this->vm().count( "help" ) ) {
        std::cout << this->optionsDescription() << "\n";
        return;
    }

    using namespace Life::vf;
    //---------------------------------------------------------------------------------------//
    
    this->changeRepository( boost::format( "ola/%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % Order
                            % this->vm()["hsize"].template as<double>()
                            );

    //---------------------------------------------------------------------------------------//

    mesh_ptrtype mesh = createMesh( meshSize );

    space_ptrtype Ah = space_type::New( mesh );
    element_type u( Ah, "u" );
    element_type v( Ah, "v" );

    element_type uexact( Ah, "uexact" );

    value_type pi = M_PI;
    AUTO( g, sin(pi*Px())*cos(pi*Py()) );

    //---------------------------------------------------------------------------------------//
    //define mesh for the curved boundaries
   
    boost::shared_ptr<struct_mesh_type> struct_mesh( new struct_mesh_type );
    GmshTensorizedDomain<1,1,1,Simplex> td;
    td.setCharacteristicLength( meshSize );
    td.setX( std::make_pair( 0, 5 ) );
    std::string fname = td.generate( "boundary_description" );
    ImporterGmsh<struct_mesh_type> struct_import( fname );
    struct_mesh->accept( struct_import );

    //---------------------------------------------------------------------------------------//
    //Define functionspace for the boundaries as well as the polynomials that define it
    
    struct_functionspace_ptrtype Xh = struct_functionspace_type::New( struct_mesh );

    struct_element_type bc_top( Xh, "bc_top" );
    AUTO( f,  val(1.0 + 0.3*cos(Px()))  );
    bc_top = vf::project( Xh, elements(Xh->mesh()), f );

    struct_element_type bc_bottom( Xh, "bc_bottom" );
    AUTO( f2, val(-1.1 - 0.3*cos(Px()))  );
    bc_bottom = vf::project( Xh, elements(Xh->mesh()), f2 );

    struct_element_type bc_reference_top( Xh, "bc_reference_top" );
    bc_reference_top = vf::project( Xh, elements(Xh->mesh()), constant(1.0) );

    struct_element_type bc_reference_bottom( Xh, "bc_reference_bottom" );
    bc_reference_bottom = vf::project( Xh, elements(Xh->mesh()), constant(-1.0) );

    struct_element_type bc_disp_top( Xh, "bc_top" );
    bc_disp_top = bc_top;
    bc_disp_top -= bc_reference_top;
    bc_disp_top.updateGlobalValues();

    struct_element_type bc_disp_bottom( Xh, "bc_bottom" );
    bc_disp_bottom = bc_bottom;
    bc_disp_bottom -= bc_reference_bottom;
    bc_disp_bottom.updateGlobalValues();

    //---------------------------------------------------------------------------------------//
    //define set of flags
    std::map< std::string, std::vector<flag_type> > flagSet;
    flagSet["fixed_bc"].push_back(3);
    flagSet["fixed_bc"].push_back(1);
    flagSet["moving_bc"].push_back(2);
    flagSet["moving_bc"].push_back(4);

    //define set of polynomials that describe boundary
    std::vector<struct_element_type> polyBoundarySet;
    polyBoundarySet.push_back(bc_disp_bottom);
    polyBoundarySet.push_back(bc_disp_top);

    std::vector<struct_element_type> referencePolyBoundarySet;
    referencePolyBoundarySet.push_back(bc_reference_bottom);
    referencePolyBoundarySet.push_back(bc_reference_top);

    //---------------------------------------------------------------------------------------//

    ALE< Simplex<2,1> > aleFactory( std::make_pair(0,5), mesh, this->vm() );

    aleFactory.generateHighOrderMap( flagSet["moving_bc"], referencePolyBoundarySet, polyBoundarySet );

    MeshMover<mesh_type> mesh_mover;
    mesh_mover.apply(mesh, aleFactory.getDisplacement() );

/*
    AUTO( fdep, vec(cst(0.),sin(pi*Px())) );
    dep_functionspace_ptrtype dXh = dep_functionspace_type::New( mesh );
    dep_element_type dep( dXh, "dep" );
    dep = project( dXh, elements(mesh), fdep );
    MeshMover<mesh_type> mesh_mover;
    mesh_mover.apply(mesh, dep );
//    uexact = project( Xh, elements(mesh), g );
//    uexact.updateGlobalValues();*/

    this->exportResults( uexact , uexact);

} // AdvReacDiff::run

template<int Dim,int Order>
void
AdvReacDiff<Dim,Order>::solve( sparse_matrix_ptrtype& D,
                       element_type& u,
                       vector_ptrtype& F )
{
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );

    M_backend->solve( D, D, U, F , false);

    u = *U;
} // AdvReacDiff::solve


template<int Dim,int Order>
void
AdvReacDiff<Dim,Order>::exportResults( element_type& U ,element_type& Uexact )
{
    Log() << "exportResults starts\n";

    exporter->step(0)->setMesh( U.functionSpace()->mesh() );

    exporter->step(0)->add( "pid",
                            regionProcess( boost::shared_ptr<p0_space_type>( new p0_space_type( U.functionSpace()->mesh() ) ) ) );
    exporter->step(0)->add( "u", U );
    exporter->step(0)->add( "uexact", Uexact );

    exporter->save();
}

// instantiation
template class AdvReacDiff<2,2>;
template class AdvReacDiff<2,5>;
